%pylab inline
import scipy
import scipy.signal
import scipy.fftpack
import numpy as np

Fs=962 # Sample rate, Hz

Populating the interactive namespace from numpy and matplotlib

RUN_NAME="foo"

import serial

s = serial.Serial("/dev/tty.usbserial-A4006Dx4", 9600)

N_SAMPLES=Fs*10

for i in range(200):
    s.read() # discard the head to discharge things

data = [ ord(s.read()) for _ in range(N_SAMPLES) ]

s.close()

np.save("%s.npy" % RUN_NAME, data)

figure(figsize=(20,5))
plot(data)
len(data)

9620

def show_fft(Fs, y):
    # sample spacing
    T = 1.0 / Fs
    # Number of samplepoints
    N = len(y)

    w = scipy.signal.blackman(N)
    ywf = scipy.fftpack.fft(y*w)

    xf = np.linspace(0.0, 1.0/(2.0*T), N/2)

    axvline(60, color="#ffcccc")
    axvline(120, color="#ccccff")


    semilogy(xf[1:N/2], 2.0/N * np.abs(ywf[1:N/2]), '-r')
    grid()
  

def show_fft_and_plot(Fs, y):
    figure(figsize=(20,5))

    subplot(1,2,1)
    show_fft(Fs,y)
    ylim(10**-3, 10**1)

    subplot(1,2,2)
    plot(y)

t0 = 4000
dt = 500

y = data[t0:t0+dt]
show_fft_and_plot(Fs, y)

lpf = scipy.signal.remez(150,[ 0.0,10,15,Fs/2], [1,0], Hz=Fs)
yp = scipy.signal.lfilter(lpf, [1], y)[len(lpf):]

figure()

#show_fft_and_plot(Fs, yp)
show_fft(Fs, yp)
xlim(0,70)
ylim(10**-5,10**-1)

(1e-05, 0.1)

hpf = scipy.signal.remez(147,[ 0.0,0.2, 0.25,0.5], [0,1])
yp = scipy.signal.lfilter(hpf, [1], y)

show_fft_and_plot(Fs, yp)

figure(figsize=(20,5))
hist(yp, bins=100)
ylim(0,200)
True

figure(figsize=(20,5))
show_fft(Fs, yp)
ylim(0.001,2)

(0.001, 2)

def gen_f(Fs, f, n):
    return np.sin(2*math.pi*f/Fs*np.linspace(0.0,n,n))


def downmix(Fs, f, y):
    return y*gen_f(Fs, f, len(y))

figure(figsize=(20,5))
#plot()

lpf_am_demod = scipy.signal.remez(50,[0,0.03,0.05,0.5], [1,0])

yd = downmix(Fs, 140, diff(y))
y_demod = scipy.signal.lfilter(lpf_am_demod, [1], abs(yd) )
#y_demod = scipy.signal.lfilter(lpf_am_demod, [1], abs(downmix(Fs, 280, yp)))  + \
#          scipy.signal.lfilter(lpf_am_demod, [1], abs(downmix(Fs, 140, yp)))  + \
#          scipy.signal.lfilter(lpf_am_demod, [1], abs(downmix(Fs, 420, yp)) )

plot(diff(y), alpha=0.2)
plot(y_demod[50:]*2, 'k')


ylim(0,300)

(0, 300)